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Abstract 

A  two-dimensional  two-phase  model  based  on  the  classical  two-fluid  model  is  used  to  analyze  electrochemical  and  thermal  transport  in  a  PEMFC. 
The  model  is  extended  to  account  for  the  dependence  of  interfacial  area  density  on  liquid  volume  fraction.  At  a  given  fixed  voltage,  the  fuel  cell 
generates  maximum  current  density  for  low  through-plane  and  high  in-plane  thermal  conductivities  at  high  humidity  operating  conditions.  It  is 
also  predicted  that  for  low  humidity  operating  conditions,  the  fuel  cell  generates  maximum  current  density  if  the  GDL  is  tailored  to  have  high 
through-plane  thermal  conductivity  near  the  inlet  and  progressively  decreasing  through-plane  thermal  conductivity  at  distances  away  from  the 
inlet.  At  fully  humidified  cathode  inlet  conditions,  narrower  current  collector  ribs  generate  higher  current  densities  at  all  voltages  by  reducing  the 
resistance  to  diffusion  of  reactants  and  products  through  the  GDL.  In  order  to  maximize  the  current  density  at  low  humidities,  ribs  must  be  wider 
near  the  inlet  and  narrower  away  from  the  inlet.  The  proposed  methodology  for  tailoring  GDL  through-plane  thermal  conductivities  and  rib  widths 
reduces  the  risk  of  membrane  dehydration  near  inlet  and  also  reduces  the  possibility  of  excessive  liquid  accumulation  in  the  region  away  form  the 
inlet. 

©  2008  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Over  the  past  few  years,  there  has  been  an  increasing  inter¬ 
est  in  the  use  of  proton  exchange  membrane  (PEM)  fuel  cells 
(FCs)  as  an  environmentally  friendly  power  source.  In  these 
cells,  oxygen  and  hydrogen  combine  to  form  water,  and  a  part 
of  the  released  energy  from  the  reaction  is  directly  converted  to 
electricity. 

Liquid  water  is  produced  in  the  cathode  catalyst  layer  due  to 
electrochemical  reactions.  Excessive  amount  of  liquid  water  in 
the  cathode  GDL  reduces  the  flow  of  reactants  to  the  catalyst 
layer  by  clogging  the  GDL  pores.  It  also  reduces  the  reaction 
rate  by  covering  the  reaction  sites  [1],  The  amount  of  liquid 
water  accumulated  in  the  cathode  GDL  depends  on  many  fac¬ 
tors,  including  the  temperature  distribution  inside  the  cell  as  it 
controls  the  phase  change  of  water.  Over  the  past  few  years  many 


*  Corresponding  author.  Tel.:  +1  814  865  1345. 

E-mail  addresses:  cjb282@psu.edu  (C.J.  Bapat),  Thynell@psu.edu 
(S.T.  Thynell). 

0378-7753/$  -  see  front  matter  ©  2008  Elsevier  B.V.  All  rights  reserved, 
doi:  10. 1016/j.jpowsour.2007. 12.033 


different  models  have  been  developed  to  analyze  the  water  and 
heat  management  problems  in  a  fuel  cell.  The  earliest  models 
[2-6]  did  not  include  an  explicit  two-phase  model  to  account 
for  the  liquid  water  transport  in  the  catalyst  layers  and  GDLs. 
Thus,  these  single-phase  models  were  unable  to  accurately  pre¬ 
dict  the  polarization  behavior  at  high  current  densities  where  the 
effects  of  liquid  water  accumulation  are  most  severe  [7].  Hence, 
different  two-phase  models  were  developed  [7-20]  to  account 
for  the  effect  of  liquid  water  transport.  The  two-phase  models 
were  based  on  the  unsaturated  flow  theory  [2 1  ]  or  the  multiphase 
mixture  model  [22],  Meng  [23]  and  Weber  and  Newman  [24,25] 
further  developed  improved  models  for  water  transport  in  the 
ionic  membrane.  The  effects  of  various  flow  and  transport  prop¬ 
erties  on  the  polarization  behavior  were  analyzed  using  these 
single-phase  and  two-phase  models. 

The  effect  of  the  thermal  conductivity  of  the  GDLs  on 
polarization  behavior  was  analyzed  by  Ju  et  al.  [26]  using  a 
single-phase  model.  The  authors  also  assumed  the  GDL  thermal 
conductivity  to  be  isotropic.  The  GDLs  are  commonly  con¬ 
structed  of  carbon  fibers  either  in  woven  or  paper  form.  Electron 
micrographs  for  different  GDLs  [27]  suggest  that  the  carbon 
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Nomenclature 

a  water  vapor  activity 

A  active  area  (m2) 

Across  area  of  cross-section  of  the  channel  (m2) 

Aig  gas-liquid  interfacial  area  per  unit  volume  (m-  1 ) 
c  specific  heat  capacity  (kJ  kg- 1  K- 1 ) 

C  molar  concentration  (mole  m-3) 

r/hyd  hydraulic  diameter  (m) 

D  diffusivity  (m2  s-1) 

E W  equivalent  weight  (kg  equiv- 1 ) 

F  Faradays  constant  (C  mole- 1 ) 

GDL  gas  diffusion  layer 

h  convective  heat  transfer  coefficient  (W  m-2  K- 1 ) 

i,  i  local  ionic  current  density  (A  m-2) 
igef  reference  current  density  (A  m-3) 

itot  total  current  density  produced  by  the  fuel  cell 

(Am-2) 

j  volumetric  electrochemical  reaction  rate  (A  m-  3 ) 

k  thermal  conductivity  ( W  m- 1  K- 1 ) 
kf  thermal  conductivity  of  fluid  (W  m-1  K-1) 
kxx  through-plane  thermal  conductivity 

(Wm-1K-1) 

kyy  in-plane  thermal  conductivity  (W  m-1  K-1) 

KgI  relative  permeability  for  gas  phase 

K\\  Henry’s  law  constant  (Pam3  mole-1) 

A-]!-  relative  permeability  for  liquid  phase 

Kv  hydrauhc  permeability  (m2) 

MEA  membrane  electrode  assembly 

noe  electro-osmotic  drag  coefficient 

Nu  Nusselt  number 

p  pressure  (Pa) 

P  perimeter  of  cross-section  (m) 

PEM  proton  exchange  membrane 

R  universal  gas  constant  (Jmole-1  K-1),  electric 
resistance  (f 2 ) 

S  volumetric  source  term  (mole  m-3  s- 1 ) 

t  time  (s) 

T  temperature  (K) 

U0  open  circuit  potential  (V) 

v  velocity  (ms-1) 

Vo  molar  specific  volume  of  water  (m3  mole-1) 

Vceii  cell  voltage  (V) 

AVnet  total  volume  of  all  the  computational  cells  adja¬ 

cent  to  the  GDL/bipolar  plate  interface 

Greek  letters 

a  phase  volume  fraction 

aa  anodic  transfer  coefficient 

ac  cathodic  transfer  coefficient 

yig  interphase  molar  transfer  coefficient  (ms-1) 

e  porosity 

9  contact  angle  (°) 

k  membrane  ionic  conductivity  (S  m-1 ) 

p.  viscosity  (Pa  s) 


p 

density  (kg  m  3) 

a 

surface  tension  (Nm-b 

),  electrical  resistivity 

(dm) 

T 

viscous  stress  (Pa) 

<t> 

potential  (V) 

Subscripts 

a 

anode 

c 

capillary,  cathode 

eq 

equivalent 

g 

gas 

1 

liquid 

m 

membrane 

P 

constant  pressure 

por 

porous  material 

wp 

wetting  phase 

Superscripts 

1 

component 

fibers  are  preferentially  oriented  in  the  in-plane  direction.  The 
thermal  conductivity  measurements  for  polymers  with  preferen¬ 
tial  fiber  orientation  [28-31]  and  also  for  certain  GDLs  [32,33] 
suggest  that  the  thermal  conductivity  of  the  GDLs  is  highly 
anisotropic.  The  effect  of  anisotropic  thermal  conductivity  on 
multicomponent  transport  in  a  fuel  cell  cathode  was  considered 
by  Pasaogullari  et  al.  [34]  using  the  multiphase  mixture  model 
[22]  and  by  Pharoah  et  al.  [35]  using  a  single-phase  model. 
In  this  paper,  we  extend  our  previous  work  [36]  to  present  an 
analysis  of  the  effect  of  both  the  magnitude  and  anisotropy  of 
the  GDL  thermal  conductivity  on  polarization  characteristics 
at  low  and  high  humidity  operating  conditions  using  a  two- 
dimensional  two-phase  model.  We  also  analyze  the  effect  of 
current  collector  rib  width  on  polarization  behavior.  We  use 
the  classical  two-fluid  model  proposed  by  Ishii  [37],  Drew  and 
Passman  [38]  and  Hassanizadeh  and  Gray  [39,40]  in  our  work 
and  also  account  for  the  effect  of  electrical  and  thermal  con¬ 
tact  resistance  at  the  GDL/bipolar  plate  interface.  In  the  next 
section  we  describe  in  detail  the  mathematical  model  used  for 
our  work. 

2.  Model  formulation 

The  model  presented  here  (see  Table  1)  accounts  for  flow 
of  reactants,  products  and  ionic  species,  as  well  as  the  genera¬ 
tion  and  transport  of  heat.  The  two-dimensional  computational 
domain  used  in  our  work  is  shown  in  Fig.  1  along  with  a 
schematic  for  a  PEMFC.  The  assumptions  and  simplifications 
used  in  our  model  are  described  next. 

2.1.  Flow  model 

We  use  a  two-phase  model  to  analyze  transport  of  species  in 
the  catalyst  layers  and  GDLs.  As  stated  in  Section  1,  the  two- 
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Table  1 

Model  equations  and  source  terms 
Equation 


Electrochemical  reactions: 

-V  •  OcV0m)  =  FSh+ 

FSK+=jc  = 

“g!ocf  [exp{^r^_  0m)} -exp{-95rW>s-</>m)}]  +“i‘ocf  [exp{^f(0s -  exp  {-^(0S  -  0m)}]  al 

cathode 
FSH+  =  h  = 

“gt'of  [exP  {lr(0s  _  0m)}  -  eXP  {_^f ' (0»  -  0m)}]  2  +“l!oaf  [exP  {"^(0*  -  0m)}  ~  exp  {-^(0*  -  0m)}] 

Gas/Liquid  phase  species  transport: 

For  t  =  02,  H2,  N2 
8mg/icL 


V  ■  (etgg/iCg/jVg/i)  =  V  ■  (g«g/i  D^VC^,)  =F  >4  (frC*  -  C{)  Alg  + 
In  catalyst  layers 

S^i  =  4/?“g/tCf  [exP  {  ^-(0i  -  0m)}  -  exp  {-^r(0s  -  0m)}] 

S™?\  =  -^“g/P'cf  [eXP  {  (0*  _  0m)}  -  exp  (0S  -  0m)}]  (cfl) 


id  =  0  everywhere  else.  Nitrogen  does  not  take  part  in  chemical  re 
Component  balance  for  water: 

In  catalyst  layers  and  GDLs 

^ ^  +  V  •  (e«g/iCg7°vg/,)  =  V  •  (eag/1D$°VC$°)  =p  y»r o  (cf*°  -  Alg  +  5^° 

In  membrane 

^a?20  =  V  '  (£“gOg2°VCg2°)  +  s7° 

Sg2°  =  — V  ■  (  yfi)  in  membrane  and  catalyst  layers  SjH2°  =  —  S’  *n  cathode  catalyst  layer 
Two  phase  momentum  balance: 

3g«s/iPg/rvs/i  +  v  .  £Q,g/]Pg/lVg/lVg/1  =  -Vsofg/ipg/i  +  V  •  eag/iTg/i  -  andpi  =pg  - pc 

Heat  balance  equation: 

|(PeqCeq7’)  +  V  •  (eamymT)  +  V  -  (eagPgygCpgf)  =  V  •  (keq  •  VT)  +  y%°  (c^°  -  AlgMn2°Lw  +  q>" 

PeqCeq  =  (1  -  £)PporCpor  +  E<X\P\Cl  +  eagPgCpg 
keq  =  (1  -  £)kpor  +  eotlkl  +  ECtgkg 

where  q'"  =  |ya(0s  —  0m)l  +  l—  in  anode  catalyst  layer  q'"  =  ~  in  membrane  q'"  =  |yc(0s  —  0m)l  +  ~  +  | 5TAS  jy  I  in  cathode  catalyst  layer 
q'"  =  tfotcr  in  GDLs  q'"  =  ('1°lA}yC’c°ll‘‘‘C‘  ^or  a  computational  cell  adjacent  to  GDL/bipolar  plate  interface  Also  see  Table  3  for  a  simplified 
expression  for  kcq. 


Ref. 


[41] 


[39,40] 


[39,40] 


[39,40] 


[36] 


phase  transport  of  different  species  is  modeled  using  the  classical 
two-fluid  model  proposed  by  Ishii  [37],  Drew  and  Passman  [38] 
and  Hassanizadeh  and  Gray  [39,40],  The  porous  drag  is  mod¬ 
eled  using  a  source  term  based  on  the  Darcy  law.  The  effects 
of  phase  change  and  species  dissolution  are  modeled  using  rate 
expressions  which  account  for  the  effect  of  gas-liquid  interfacial 
area  density. 

We  use  a  single-phase  model  for  transport  inside  the  mem¬ 
brane  which  means  that  we  do  not  allow  the  liquid  water 
produced  at  cathode  to  enter  the  membrane  and  consider  only  the 
concentration-driven  diffusion  of  water  vapor  through  the  mem¬ 
brane.  We  also  do  not  consider  phase  change  of  water  inside  the 
membrane.  The  assumption  of  single-phase  transport  inside  the 
ionic  membrane  seems  reasonable  for  two  reasons:  the  hydraulic 
permeability  of  the  membrane  is  a  few  orders  of  magnitude  lower 
than  the  hydraulic  permeability  of  the  GDL  (see  Table  2)  and  the 
membrane  is  known  to  be  at  least  slightly  hydrophobic  with  pore 
radii  much  smaller  than  the  GDL  [24,25],  The  overall  effect  of 


Table  2 

Model  parameters  and  properties 


C.J.  Bapat,  S.T.  Thynell  /Journal  of  Power  Sources  179  (2008)  240-251 


243 


Quantity 

Value 

Gas  channel  length  (m) 

1.577  [54] 

Gas  channel  width  (m) 

2.16E-3  [54] 

Current  collector  width  (m) 

0.89E-3  [54] 

Active  area  (nr) 

50E-4  [54] 

Membrane 

Nafion®  112  [54] 

Gas-diffusion  layer  thickness  (m) 

350E-6  [48] 

Catalyst  layer  thickness  (m) 

10.0E-6 

Height  of  the  symmetrical  section  (m) 

1.525E— 3 

Anode  humidification  temperature  (K) 

353  [54] 

Cathode  humidification  temperature  (K) 

353  [54] 

Anode  inlet  stoichiometry  (A  m-2  equiv) 

1.875E4  [54] 

Cathode  inlet  stoichiometry  (A  m-2  equiv) 

1.125E4  [54] 

Coolant  water  temperature  (K) 

353  (assumed) 

Cathode  inlet  pressure  (atm) 

1.5  [54] 

Anode  inlet  pressure  (atm) 

1.5  [54] 

Gas-diffusion  layer  electrical  resistivity  (£2  m) 

80.0E-5  [32] 

Through-plane  thermal  conductivity  of  gas-diffusion  layers  k **  (W  m-1  K-1) 

0.22  [49] 

Ratio  of  thermal  conductivities  kyy/kx* 

20:1  (assumed  based  on  [32]) 

Thermal  conductivity  of  the  membrane  (W  m-1  K-1) 

0.16  [49] 

Thermal  conductivity  of  the  catalyst  layer  (W  m-1  K-1) 

0.27  [49] 

Thermal  conductivity  of  current  collector,  steel  (W  m-1  K-1) 

16.0 

GDL/Bipolar  plate  thermal  contact  conductance  (W  m-2  K-1 ) 

10,000  (assumed  based  on  [61]) 

GDL/Bipolar  plate  electrical  contact  resistance  (mt!  cm1 2) 

41.7  [54] 

Gas-diffusion  layer  porosity 

0.74  [27] 

Maximum  interfacial  area  density  (m-1) 

500 

Catalyst  layer  porosity 

0.15  (fitted) 

Hydraulic  permeability  of  GDL  (m2) 

6.3E— 12  [48] 

Hydraulic  permeability  of  membrane  (m2) 

1.0E-18  [25] 

Hydraulic  permeability  of  catalyst  layer  (m2) 

1.0E— 13  (assumed  based  on  [48]) 

Contact  angle  in  GDL  and  catalyst  layer  (°) 

120 

Hydrogen  diffusivity  (m2  s-1) 

0.915E— 4  [62] 

Oxygen  diffusivity  (m2  s-1) 

3.0E-5  [62] 

Water  vapor  diffusivity  (m2  s-1) 

3.0E-5  [62] 

Electro-osmotic  drag  coefficient 

1.0  [63] 

Transfer  coefficient  at  anode,  a a  =  ac 

1.0  [20] 

Anode  reference  exchange  current  density  (Am-3) 

1.5E9  (assumed) 

Transfer  coefficient  at  cathode,  a a  =  ac 

1.0  (Vceii<0.5V)[51] 

0.5  (VCeii>0.5V) 

Cathode  reference  exchange  current  density  (Am-3) 

1.1E2  (Vceii  <0.5  V)  (fitted) 

3.0E5  (Vceu>0.5V) 

Entropy  change  AS  for  Ho  +  ]  O2  — *■  H20(iiq)  (J  mol-1  K-1) 

-162.4 

using  a  single-phase  model  in  the  membrane  is  to  slightly  overes¬ 
timate  the  effects  of  liquid  water  accumulation  in  catalyst  layers 
and  GDLs  on  polarization  behavior. 

2.2.  Heat  transport  model 

We  also  assume  that  the  gas  and  liquid  phases  as  well  as  the 
porous  matrix  are  in  local  thermal  equilibrium.  Hence  we  use 
a  single  equation  derived  by  summation  of  the  respective  phase 
heat  transport  equations  to  estimate  the  temperature  distribution. 
The  heat  generation  is  modeled  using  appropriate  source  terms 
in  different  regions  of  the  fuel  cell.  The  model  accounts  for  heat 
generation  due  to 

1.  activation  polarization  in  catalyst  layers, 

2.  ohmic  polarization  in  catalyst  layers  and  the  ionic  membrane, 


3.  reversible  heat  loss, 

4.  phase  change  and 

5.  heat  generation  due  to  flow  of  electric  current  in  the  GDLs 
and  at  the  GDL/bipolar  plate  interface. 

Also,  the  entire  reversible  heat  loss  is  assumed  to  occur  in  the 
cathode  catalyst  layer. 

2.3.  Electrochemistry 

We  solve  for  transport  of  protons  in  the  catalyst  layers  and 
membrane.  The  solid-phase  potential  in  the  catalyst  layers  is 
assigned  according  to  the  boundary  conditions.  The  effect  of 
electrical  contact  resistance  is  also  taken  into  account  in  the 
boundary  conditions.  The  electrochemical  reactions  are  modeled 
by  the  Butler- Volmer  reaction  rate  expressions  [41]  which  are 
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modified  to  include  the  effect  of  liquid  accumulation  on  the 
reaction  rates  by  assuming  proportionality  to  the  phase  volume 
fractions. 

2.4.  Channel  flow  model 

We  solve  for  pressures  and  component  balances  in  the  chan¬ 
nels  to  determine  the  pressures  and  concentrations  of  different 
species  at  different  distances  form  the  inlet.  The  pressures  and 
concentrations  so  determined  are  then  used  as  boundary  con¬ 
ditions  for  the  two-dimensional  model  used  in  our  work.  This 
allows  us  to  simulate  the  polarization  behavior  at  different  dis¬ 
tances  form  the  inlet. 

The  pressure  variation  is  calculated  using  the  Darcy  friction 
factor  law  [42]  to  account  for  the  effect  of  wall  friction.  The 
reactant  concentrations  are  calculated  using  simple  mass  bal¬ 
ances  accounting  for  the  production/consumption  of  reactants 
and  products  [43].  The  water  produced  at  the  cathode  is  assumed 
to  mix  with  the  cathode  gas  stream  in  the  form  of  water  vapor, 
and  phase  change  of  water  vapor  in  the  gas  channels  is  neglected. 
The  temperature  change  of  the  reactant  gases  is  also  neglected. 
In  spite  of  these  simplifications,  the  above  method  of  calculat¬ 
ing  the  parameters  takes  into  account  the  effect  of  reactant  flow 
rates,  channel  length  and  active  area,  which  otherwise  cannot  be 
accounted  for  in  a  two-dimensional  model.  It  is  noticed  that  the 
parameter  values  calculated  above  depend  on  the  current  den¬ 
sity,  which  is  obtained  from  the  solution  of  the  electrochemistry 
model.  Hence  the  pressure  and  concentration  boundary  condi¬ 
tions  are  also  updated  during  each  iteration  of  the  numerical 
scheme. 

3.  Constitutive  relations 


4.  Boundary  conditions 

The  boundary  conditions  required  at  various  interfaces  are 
now  described  in  detail. 


4.1.  Flow  and  concentration  boundary  conditions 

The  channel  pressures  as  well  as  the  concentrations  of  reac¬ 
tants  and  products  on  the  anode  and  cathode  sides  are  specified 
as  boundary  conditions  at  the  GDL/gas  channel  interfaces.  Pres¬ 
sures  and  concentrations  at  different  distances  from  the  inlet  are 
calculated  using  the  channel  flow  model  described  previously. 
As  we  use  a  single-phase  model  in  the  ionic  membrane,  the  liq¬ 
uid  velocity  at  the  catalyst  layer/membrane  interface  is  set  to 
zero. 


In  order  to  complete  the  description  of  the  model  presented 
in  the  previous  section,  constitutive  relations  for  various  proper¬ 
ties  are  required.  The  constitutive  relations  used  in  this  work  are 
given  in  Table  3.  We  developed  a  correlation  for  the  interfacial 
area  density  using  the  existing  pore-scale  modeling  and  experi¬ 
mental  data  [44—46]  for  different  porous  media.  The  correlation 
is  plotted  in  Fig.  2  along  with  the  experimental  and  modeling 
data,  and  it  shows  an  intermediate  shape  between  the  model 
results.  We  also  correlated  the  available  experimental  data  [47] 
for  diffusivity  of  water  vapor  through  the  ionic  membrane.  We 
use  the  correlation  developed  by  Gostick  et  al.  [48]  for  capillary 
pressure  in  a  GDL.  The  E-Tek  carbon  cloth  used  by  Gostick 
and  coworkers  for  their  measurements  was  almost  completely 
hydrophilic,  whereas  experiments  show  that  the  carbon  cloth 
used  in  the  E-Tek  electrode  (E-Tek  electrode  has  been  used  in 
the  experimental  study  used  to  validate  the  model)  is  hydropho¬ 
bic  [27].  Therefore  we  correct  the  capillary  pressure  obtained  by 
their  correlation  for  contact  angle.  We  also  express  the  equiva¬ 
lent  thermal  conductivity  of  the  GDL  (see  Table  1)  in  terms  of 
the  thermal  conductivity  of  a  dry  GDL  to  facilitate  the  use  of 
existing  experimental  data  [49].  We  use  published  correlations 
for  membrane  water  content  at  80 °C  [50],  ionic  conductiv¬ 
ity  [25]  and  temperature  dependence  of  cathode  reaction  rate 
[51]. 


4.2.  Electrochemistry 

In  the  current  model  the  electrochemical  equations  are  solved 
over  the  MEA.  At  the  GDL/anode  interface,  they  are 


fls  —  U0  —  Lee]]  —  JtotA/?a, contact  —  fiot^  /?gd]  (la) 

V0m  =  0  (lb) 

At  the  GDL/cathode  interface,  they  are 

</>s  =  *totA/?c,  contact  +  hot  4  A’gdl  (2a) 

V0m  =  0  (2b) 


Again  the  analysis  presented  above  is  simplified  and  more  com¬ 
plicated  issues  involving  the  anisotropic  electrical  conductivity 
of  the  GDLs  [32,33]  have  not  been  addressed.  It  should  also  be 
noted  that  the  above  boundary  conditions  are  dependent  on  the 
current  density.  They  are  also  updated  during  each  iteration  of 
the  numerical  scheme,  as  in  case  of  the  flow  boundary  condi¬ 
tions.  The  corresponding  heat  generation  terms  due  to  contact 
and  GDL  electrical  resistance  are  also  added  to  the  heat  equation 
in  the  GDLs. 
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Table  3 

Constitutive  relations 


Expression 

Ref. 

Capillary  pressure: 

“wr  =  (l+(£fr  , 

where/=  (^)  Z  and/c  =0.6982,  n  =  3.465  and  m  =  0.7114;pc  =_pc,aii-water  cos0 

[48] 

Corrections  in  property  values: 

D‘g  =  Og  0e1'5Q!g0;  Kg,  =  a30,  Kb  =  af  ° 

Interfacial  area  density: 

A,g  =  2MAlgm^a°x6a\2 

See  Fig.  2 

Membrane  water  content  at  80  °C: 

X  =  0.300  +  10.8a  -  16.0a2  +  14.1a3 

[50] 

Ionic  conductivity: 

K  =  0.5 (/  -  0.06)15  exp  -  f )] 

where  /  =  y^-uvo  aRd  Vm  =  ^ 

[25] 

Vapor  diffusivity  in  the  membrane: 

^2mem  =  (3.18633  x  10“10)a  -  5.23356  x  10-11 

Temperature  dependence  of  exchange  current  density: 

‘ocf(D  =  i'S%(353K)exp  [-TT  (t  -  353T5)] 

[51] 

Equivalent  thermal  conductivity: 

The  expression  for  keq  given  in  Table  1  can  be  simplified  as,  keq  =  kjry  +  Ecufki  —  kg)  where  k^ry  =  |  q”  ^  j  and  k and  kyy 

[49] 

are  experimentally  measured  through-plane  and  in-plane  thermal  conductivities.  The  liquid  and  gas  thermal  conductivities  are 
isotropic. 

4.3.  Heat  transfer 

The  heat  transfer  at  the  GDL/gas  channel  interface,  current 
collector/gas  channel  interface  and  the  current  collector/coolant 
channel  interface  occurs  due  to  convection.  Characteristic  Re 
numbers  for  flow  in  the  gas  and  coolant  channels  are  low  and 
in  the  laminar  range.  Further,  the  hydrodynamic  and  thermal 
entry  lengths  are  relatively  short  compared  to  channel  lengths. 
Hence,  the  convective  heat  transfer  coefficients  for  the  gas  and 
coolant  channels  are  calculated  using  the  following  relation 
applicable  for  laminar  fully  developed  flows  in  square  channels 
[52], 

Nu=hdb^  =  2.98  (3.1) 

kf 

where 

4  a 

4  yd  =  — (3.2) 

Additionally,  the  channel  boundaries  are  at  nearly  a  uniform 
temperature.  The  effect  of  a  variable  temperature  along  the  GDL 
boundary  on  the  overall  heat  transfer  coefficient  is  neglected. 
The  gas  and  coolant  temperatures  were  assumed  to  be  at  353  K 
(cell  operating  temperature). 

Also,  as  the  two-dimensional  domain  used  in  our  work  is  a 
symmetrical  section,  all  other  boundaries  are  assumed  to  obey 
symmetry  boundary  conditions. 


5.  Numerical  scheme 

The  set  of  two-dimensional  conservation  equations  for 
species,  electrochemistry,  fluid  flow  and  heat  conduction  was 
solved  iteratively  over  the  computational  domain  shown  by  dot¬ 
ted  lines  in  Fig.  1.  The  values  of  various  parameters  used  for 
the  model  are  given  in  Table  2.  The  heat  conduction  equation 
was  solved  until  the  relative  error  in  the  overall  heat  balance 
was  <10-3,  and  the  flow  equations  were  solved  using  the  SIM¬ 
PLER  scheme  [53]  until  the  relative  error  mass  balance  was 
<10_1°.  The  iterations  continued  until  convergence  criteria  for 
all  the  equations  were  simultaneously  met.  The  solution  was 
obtained  for  successively  refined  grids  and  a  Cartesian  grid  of 
size  31  x  10  provided  a  grid  independent  solution.  The  polar¬ 
ization  curve  was  obtained  using  the  model,  and  the  cathode 
reference  exchange  current  density  was  used  as  a  fitting  param¬ 
eter  [4 1  ]  along  with  the  catalyst  layer  porosity.  Since  the  porosity 
of  the  E-Tek  electrode  catalyst  layer  used  in  the  experiment  [54] 
used  for  validation  is  not  known,  the  electrode  porosity  was  used 
as  a  fitting  parameter  for  the  model.  Mench  et  al.  [54]  humidify 
the  anode  stream  to  a  temperature  above  the  fuel  cell  operating 
temperature.  Since  the  fuel  cell  is  maintained  at  a  lower  tem¬ 
perature  than  the  anode  humidification  temperature,  the  excess 
water  vapor  in  the  anode  stream  will  most  probably  condense 
near  the  inlet.  It  is  not  possible  to  predict  using  our  model  if 
the  water  vapor  condenses  and  is  carried  forward  in  the  form  of 
droplets  in  the  gas  stream  or  if  the  water  vapor  condenses  at  the 
anode  GDL/gas  channel  interface  and  then  penetrates  the  anode 
GDL.  Hence  for  the  purposes  of  this  model  we  assume  the  anode 
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to  be  fully  humidified  at  the  fuel  cell  temperature  and  that  none 
of  the  condensed  water  penetrates  the  anode  GDL. 

The  experimental  polarization  curves  at  different  distances 
from  the  inlet  [54]  along  with  the  model  predictions  are  shown 
in  Fig.  3.  It  can  be  seen  that  even  with  a  two-dimensional  model 
the  polarization  behavior  at  different  distances  from  the  inlet  can 
be  predicted  with  a  good  degree  of  accuracy. 

6.  Discussion  of  results 

In  this  work,  we  study  the  effects  of  the  anisotropic  ther¬ 
mal  conductivity  of  the  GDLs,  GDL  width  and  the  width  of 
the  current  collector  ribs  on  the  polarization  behavior  at  differ¬ 
ent  operating  conditions.  Fig.  3  shows  that  the  effect  of  liquid 
accumulation  in  the  cathode  GDL  and  catalyst  layer  is  more 
significant  at  lower  voltages.  In  general,  the  operating  voltage 
depends  on  the  type  of  application  for  which  the  fuel  cell  is 
used.  Higher  voltages  generate  lower  current  densities  at  higher 
efficiencies,  whereas  lower  voltages  generate  higher  current 
densities  but  at  lower  efficiencies.  We  analyze  results  at  three 
different  voltages  (0.65,  0.6  and  0.55  V)  at  which  the  current 
densities  range  from  about  50-80%  of  the  maximum  current 
density  in  Fig.  3.  Also,  the  inlet  flows  are  generally  humidi¬ 
fied  to  maintain  the  ionic  conductivity  of  the  ionic  membrane. 
Therefore,  we  also  investigate  the  polarization  behavior  at  dif¬ 
ferent  cathode  inlet  humidities.  Also,  for  Figs.  4-10,  values  of 
the  parameters  kept  constant  are  given  in  Table  2. 

6.1.  Fully  humidified  cathode  inlet 

Figs.  4  and  5  show  the  variations  of  current  density,  max¬ 
imum  fuel  cell  temperature,  membrane  ionic  conductivity  and 
cathode  catalyst  layer  liquid  volume  fraction  near  the  inlet  for 
different  values  of  in-plane  and  through-plane  thermal  conduc¬ 
tivities  of  the  GDL.  Lower  through-plane  thermal  conductivity 
leads  to  higher  temperatures.  Higher  temperatures  allow  more 
liquid  water  to  evaporate  which  reduces  the  amount  of  liquid 
water  accumulated  in  the  cathode  catalyst  layer.  Higher  tem¬ 
peratures,  however,  also  reduce  the  water  vapor  activity  and 


Current  Density  (A  m  2) 

Fig.  3.  Comparison  of  experimental  [54]  and  model  polarization  results. 


Fig.  4.  Variation  of  (a)  current  density  and  (b)  ionic  conaucuvity  with  thermal 
conductivity  of  the  GDL  near  the  inlet  for  fully  humidified  cathode  inlet  flow. 


hence  adversely  affect  the  membrane  ionic  conductivity.  Thus 
increasing  the  through-plane  thermal  conductivity  of  the  GDL 
leads  to  two  competing  effects:  an  increase  in  the  amount  of 
liquid  water  at  the  cathode  catalyst  layer  tends  to  reduce  the 
current  density,  whereas  a  higher  ionic  conductivity  of  the 
membrane  tends  to  increase  the  current  density.  Since  the  rela¬ 
tive  magnitudes  of  these  two  competing  effects  depend  on  the 
through-plane  thermal  conductivity  of  the  GDL,  the  variation  of 
current  density  with  through-plane  thermal  conductivity  can  be 
expected  to  show  a  maximum  at  a  particular  value  of  the  through- 
plane  thermal  conductivity.  The  variation  of  current  density 
with  increasing  through-plane  thermal  conductivity  (see  Fig.  4a) 
shows  a  maximum  for  low  in-plane  thermal  conductivities  at 
lower  voltages.  At  higher  in-plane  thermal  conductivities,  how¬ 
ever,  the  current  density  decreases  with  increasing  through-plane 
thermal  conductivity  for  all  the  voltages.  Thus,  at  higher  in-plane 
thermal  conductivities,  the  adverse  effect  of  higher  liquid  accu¬ 
mulation  in  the  catalyst  layer  dominates  and  the  current  density 
decreases  with  increasing  through-plane  thermal  conductivity. 
At  low  in-plane  thermal  conductivities  and  lower  voltages  the 
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Fig.  5.  Variation  of  maximum  fuel  cell  temperature  and  catalyst  layer  liquid 
volume  fraction  with  thermal  conductivity  of  the  GDL  near  the  inlet  for  fully 
humidified  cathode  inlet  flow. 

temperatures  are  higher  and  the  membrane  ionic  conductivities 
are  lower.  Hence  the  variation  of  current  density  with  increas¬ 
ing  through-plane  thermal  conductivity  shows  a  maximum  at 
lower  voltages  and  lower  in-plane  thermal  conductivities.  Lower 
voltages  generate  higher  current  densities  which  lead  to  the  pro¬ 
duction  of  more  liquid  water  and  also  higher  temperatures  due 
to  higher  heat  generation.  Hence  it  is  important  to  determine  if 
the  higher  temperatures  produced  at  lower  voltages  are  able  to 


(a)  (b) 


Fig.  6.  Comparison  of  temperature  distributions  in  the  cathode  GDL  for  low 
and  high  in-plane  GDL  thermal  conductivities  for  fully  humidified  cathode  inlet 
flow,  (a)  kxx  =  l.OWm-1  K-1  and  kyy/kxx  =  1.0;  (b)  kxx  =  l.OWm"1  K"1  and 

kyy/kxx  =  10. 
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Cathode  Humidification  Temperature  (K) 


Fig.  7.  Variation  of  current  density  with  cathode  humidification  temperature  at 
different  distances  from  inlet. 


320  325  330  335  340  345  350  355 
Cathode  Humidification  Temperature  (K) 


Fig.  8.  Variation  of  (a)  catalyst  layer  liquid  volume  fraction  and  (b)  ionic  con¬ 
ductivity  with  cathode  humidification  temperature  at  different  distances  from 
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Cathode  Humidification  Temperature  (K) 


Fig.  9 .  Variation  of  current  density  with  cathode  humidification  temperature  near 
fuel  cell  inlet  for  different  through-plane  thermal  conductivities  of  the  GDL. 


lqx  (Wnf'lC1) 


plane  thermal  conductivity  of  the  GDL  for  different  rib  widths  near  inlet  at  fully 
humidified  cathode  inlet  flow  condition. 


evaporate  the  higher  amount  of  liquid  water  produced  at  lower 
voltages.  Fig.  5  shows  that,  at  lower  through-plane  thermal  con¬ 
ductivities  the  predicted  maximum  temperatures  are  higher  for 
lower  voltages.  Because  of  these  high  temperatures,  the  liquid 
water  is  able  to  evaporate  and  liquid  accumulation  in  the  cath¬ 
ode  catalyst  layer  is  less  for  lower  voltages  than  that  for  higher 
voltages.  However,  as  the  through-plane  thermal  conductivity 
increases,  the  temperatures  are  lowered  and  the  liquid  accumu¬ 
lation  in  the  cathode  catalyst  layer  at  lower  voltages  increases 
to  a  value  above  that  for  higher  voltages.  Even  though  higher 
temperatures  reduce  liquid  accumulation  in  the  catalyst  layer 
and  GDL,  such  high  temperatures  also  tend  to  adversely  affect 
the  long  term  durability  of  the  ionic  membrane  [55],  Hence  it 
is  desirable  to  achieve  higher  current  densities  without  creating 
excessively  high  temperatures. 

The  effects  of  increasing  in-plane  thermal  conductivity  on 
current  density,  maximum  fuel  cell  temperature,  membrane  ionic 
conductivity  and  cathode  catalyst  layer  liquid  volume  fraction 
are  also  shown  in  Figs.  4  and  5.  Higher  in-plane  thermal  con¬ 
ductivities  lead  to  lower  temperatures.  Lower  temperatures  lead 
to  higher  liquid  volume  fractions  in  the  cathode  catalyst  layer 
and  also  higher  ionic  conductivity.  It  can  also  be  noticed  that 
the  current  density  remains  almost  the  same  for  low  through- 
plane  thermal  conductivities  with  increasing  in-plane  thermal 
conductivity  for  all  voltages  shown  in  Fig.  4a.  As  shown  in  Fig.  6, 
higher  in-plane  thermal  conductivities  create  more  uniform  tem¬ 
peratures  than  lower  in-plane  thermal  conductivities.  It  can  also 
be  noticed  in  Fig.  6  that  the  maximum  temperature  is  lower  if 
the  in-plane  thermal  conductivity  of  the  GDL  is  high.  Higher 
in-plane  thermal  conductivities  thus  reduce  the  risk  of  forma¬ 
tion  of  hot  spots  in  the  region  directly  exposed  to  gas  channel 
and  also  reduce  the  possibility  of  excessively  low  temperatures 
in  the  region  directly  in  contact  with  the  current  collector.  More 
uniform  temperatures  thus  reduce  the  susceptibility  of  local  dry¬ 
out  of  the  ionic  membrane  in  the  region  directly  exposed  to  gas 
channel  and  also  reduce  the  possibility  of  excessive  liquid  accu¬ 
mulation  in  the  low  temperature  regions.  Thus  a  combination 
of  high  in-plane  thermal  conductivity  and  a  low  through-plane 
thermal  conductivity  of  the  GDL  yields  higher  current  densities. 

The  effects  of  liquid  water  accumulation  are  more  severe  in 
regions  further  away  from  the  inlet  due  to  reduction  in  reactant 
concentration  and  a  highly  humidified  cathode  stream  created 
by  the  product  water  joining  the  cathode  gas  channel.  Hence  in 
general,  in  order  to  maximize  the  current  density  under  fully 
humidified  inlet  conditions,  the  through-plane  thermal  conduc¬ 
tivity  of  the  GDL  should  be  as  low  as  possible  and  the  in-plane 
thermal  conductivity  should  be  as  high  as  possible. 

6.2.  Low  humidity  cathode  inlet 

The  effects  of  low  humidity  operation  on  current  density, 
ionic  conductivity  and  catalyst  layer  liquid  volume  fraction  are 
shown  in  Figs.  7  and  8.  We  show  the  results  at  two  differ¬ 
ent  distances  from  the  inlet.  The  ionic  conductivity  and  liquid 
accumulation  in  the  catalyst  layer  are  reduced  due  to  the  lower 
cathode  humidity.  However,  the  reduction  in  liquid  accumula¬ 
tion  in  the  cathode  catalyst  layer  has  a  significant  impact  on 
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current  density  and  hence  the  current  density  increases  with 
lower  humidity.  If  the  inlet  cathode  humidity  is  reduced  below 
a  certain  value,  the  effect  of  lower  ionic  conductivity  starts 
to  dominate  and  the  current  density  decreases.  Thus,  for  low 
humidity  operation,  the  region  near  the  inlet  may  experience 
membrane  dehydration  while  the  regions  away  from  the  inlet 
generate  higher  current  densities.  Thus,  even  though  a  low  cath¬ 
ode  humidity  improves  current  densities  away  from  the  inlet, 
the  adverse  effect  of  membrane  dehydration  prevents  the  use  of 
low  cathode  inlet  humidities. 

In  order  to  address  the  problem  of  membrane  dehydration 
near  inlet  at  low  humidities,  Dong  et  al.  [56]  suggest  using  a 
fully  humidified  anode  inlet  flow  at  a  temperature  above  the  fuel 
cell  operating  temperature.  With  an  aim  to  develop  an  alternative 
methodology  to  address  the  problem  of  membrane  dehydra¬ 
tion  near  inlet,  we  now  investigate  the  effect  of  GDL  thermal 
conductivities  on  low  humidity  operation.  Fig.  9  shows  the  vari¬ 
ation  of  current  density  near  the  inlet  for  different  through-plane 
thermal  conductivities  at  low  humidities.  Higher  through-plane 
thermal  conductivities  result  in  higher  current  densities  at  lower 
cathode  humidities.  Higher  through-plane  thermal  conductiv¬ 
ities  in  conjunction  with  high  in-plane  thermal  conductivities 
(see  Table  2)  result  in  lower  and  more  uniform  temperatures 
inside  the  fuel  cell  which  in  turn  lead  to  a  higher  water  vapor 
activity  and  hence  reduce  the  risk  of  membrane  dehydration 
near  the  inlet.  At  lower  voltages,  higher  through-plane  thermal 
conductivity  leads  to  higher  liquid  accumulation  in  the  cath¬ 
ode  catalyst  layer.  Therefore,  the  maximum  current  density  at 
a  given  through-plane  thermal  conductivity  occurs  at  a  lower 
humidification  temperature  for  lower  voltages. 

The  analysis  presented  in  Section  6.1  shows  that  high 
through-plane  thermal  conductivities  promote  liquid  accumula¬ 
tion  in  the  catalyst  layer  at  fully  humidified  inlet  conditions.  The 
effect  of  liquid  accumulation  is  more  severe  in  regions  away  from 
the  inlet  because  the  cathode  gas  channel  humidity  increases 
towards  the  flow  channel  exit.  Therefore,  the  current  density  can 
be  increased  if  membrane  dehydration  can  be  prevented  near  the 
inlet  without  allowing  liquid  water  to  accumulate  in  the  region 
away  from  the  inlet.  The  current  density  can  thus  be  increased 
by  constructing  a  GDL  with  a  high  through-plane  thermal  con¬ 
ductivity  near  the  inlet  of  the  flow  channel  and  progressively, 
lower  through-plane  thermal  conductivity  in  regions  away  from 
the  inlet.  The  through-plane  thermal  conductivity  of  the  GDL 
at  any  distance  from  the  inlet  can  be  optimized  to  provide  the 
maximum  current  density  at  a  given  cathode  humidity  and  oper¬ 
ating  voltage.  The  GDL  can  then  be  constructed  to  have  these 
optimized  thermal  conductivities  at  appropriate  distances  along 
the  flow  channel  to  maximize  the  current  density  at  the  oper¬ 
ating  voltage.  The  through-plane  thermal  conductivity  of  the 
GDL  can  thus  be  tailored  to  reduce  the  susceptibility  of  mem¬ 
brane  dehydration  near  inlet  and  also  to  reduce  the  risk  of  liquid 
accumulation  in  the  catalyst  layer  in  regions  away  from  the  inlet. 

6,3.  Effect  of  rib  width 

The  effect  of  rib  width  variation  on  fuel  cell  operation  has 
been  analyzed  previously  [57-60]  from  the  point  of  view  of 


uid  volume  fraction  with  through-plane  thermal  conductivity  of  the  GDL  for 
different  rib  widths  near  the  inlet  at  fully  humidified  cathode  inlet  flow  condition. 

reducing  the  resistance  to  reactant  diffusion  and  electron  trans¬ 
port.  In  our  previous  work  [36]  we  used  a  single-phase  model 
to  analyze  the  effect  of  variation  of  rib  width  on  tempera¬ 
ture  distribution  inside  a  fuel  cell.  We  showed  that  narrower 
ribs  lead  to  higher  temperatures  and  using  GDLs  with  highly 
anisotropic  thermal  conductivity  (higher  in-plane  thermal  con¬ 
ductivity)  leads  to  uniform  and  lower  temperatures.  We  extend 
our  analysis  further  in  this  work. 

Figs.  10  and  1 1  show  the  effect  of  variation  of  through-plane 
thermal  conductivity  on  current  density,  ionic  conductivity,  max¬ 
imum  fuel  cell  temperature  and  cathode  catalyst  layer  liquid 
accumulation  for  different  rib  widths.  The  results  shown  in 
Figs.  10  and  1 1  are  for  fully  humidified  cathode  and  anode  flows 
near  the  inlet.  Narrower  ribs  reduce  the  resistance  to  diffusion  of 
water  vapor  from  cathode  catalyst  layer  to  the  gas  channel.  This 
reduction  in  water  vapor  diffusion  resistance  lowers  the  accu¬ 
mulation  of  liquid  water  in  the  cathode  catalyst  layer.  Hence,  for 
the  voltages  shown  in  Figs.  10  and  11,  smaller  rib  widths  lead  to 
higher  current  densities.  Fig.  1 1  also  shows  that  the  temperatures 
are  higher  for  narrower  ribs  for  all  the  considered  voltages.  Nar¬ 
rower  ribs  lead  to  higher  temperatures,  because  higher  current 
densities  lead  to  more  heat  generation  and  smaller  rib  widths 
reduce  the  area  of  heat  transfer  from  the  GDL  to  the  coolant 
channel.  Also,  the  variation  of  current  density  with  increas¬ 
ing  through-plane  thermal  conductivity  shows  a  maximum  in 
case  of  small  rib  widths  (rib : channel  =  1:5  in  Fig.  10).  Lower 
through-plane  thermal  conductivities  lead  to  higher  tempera¬ 
tures  which  reduce  liquid  accumulation  in  the  cathode  catalyst 
layer  and  also  reduce  the  ionic  conductivity  of  the  membrane. 
In  addition,  as  seen  from  Fig.  11,  the  accumulation  of  liquid 
water  in  the  cathode  catalyst  layer  is  reduced  for  narrower  ribs. 
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Fig.  12.  Variation  of  current  density  with  cathode  humidification  temperature 
near  fuel  cell  inlet  for  different  rib  widths. 


In  case  of  fuel  cells  using  narrow  ribs  and  GDLs  with  low 
through-plane  thermal  conductivities,  the  reduced  membrane 
ionic  conductivity  has  a  larger  negative  impact  on  current  den¬ 
sity  than  the  positive  effect  of  reduction  in  cathode  catalyst 
layer  liquid  accumulation.  Therefore,  for  fuel  cells  using  nar¬ 
row  ribs  and  operating  at  fully  humidified  inlet  condition,  the 
through-plane  thermal  conductivity  should  be  maintained  at  a 
higher  value  in  order  to  maximize  the  current  density  near  the 
inlet.  In  regions  further  away  form  the  inlet,  the  humidity  of  the 
cathode  gas  channel  flow  is  higher  and  lower  through-plane  ther¬ 
mal  conductivities  can  be  used  without  adversely  affecting  the 
current  density.  Therefore,  the  methodology  to  tailor  the  GDL 
through-plane  thermal  conductivity  proposed  in  Section  6.2  can 
also  be  used  to  maximize  the  current  density  for  fuel  cells  using 
narrow  ribs  and  operating  at  fully  humidified  inlet  conditions. 
In  the  present  analysis,  we  have  not  considered  the  effect  of 
reducing  the  rib  width  on  electron  transport  [58]  and  also  on 
the  structural  strength  of  the  current  collector  plate.  The  cur¬ 
rent  density  of  the  fuel  cell  operating  at  fully  humidified  inlet 
conditions  is  maximum  if  the  current  collector  rib  is  as  narrow 
as  permitted  by  the  structural  strength  and  electron  transport 
limitations. 

Fig.  12  shows  the  effect  of  rib  width  on  low  humidity 
operation  near  the  inlet.  Narrower  ribs  generate  higher  cur¬ 
rent  densities  even  at  lower  cathode  humidities  by  reducing 
the  resistance  to  diffusion  of  gases  through  the  GDL.  How¬ 
ever,  at  cathode  humidities  lower  than  an  optimum  value,  wider 
ribs  are  able  to  generate  higher  current  densities.  Wider  ribs 
are  therefore  more  effective  in  reducing  the  susceptibility  of 
membrane  dehydration  near  the  inlet.  However,  as  discussed 
earlier  in  this  section,  narrower  ribs  are  able  to  generate  higher 
current  densities  at  fully  humidified  inlet  conditions.  Hence  it 
would  be  desirable  to  use  wider  ribs  near  the  inlet  and  narrower 
ribs  in  the  region  away  from  the  inlet  in  order  to  maximize 
the  current  density  at  low  humidity  operating  conditions.  The 
width  of  the  ribs  can  thus  be  tailored  to  obtain  maximum  cur¬ 
rent  density  at  a  given  voltage  and  given  cathode  inlet  humi¬ 
dity. 


6.4.  Effect  of  GDL  thickness 

The  effect  of  reducing  GDL  thickness  on  current  density  is 
similar  to  the  effect  of  narrower  current  collector  ribs  on  current 
density.  Within  the  computational  domain  of  Fig.  1,  a  certain 
fraction  of  the  GDL  surface  is  exposed  to  the  gas  channel  and 
the  other  fraction  is  directly  in  contact  with  the  current  col¬ 
lector  ribs.  In  the  region  of  GDL  directly  exposed  to  the  gas 
channel,  smaller  GDL  thickness  reduces  the  resistance  to  diffu¬ 
sion  of  water  vapor  and  reactants  through  the  GDL.  However, 
smaller  GDL  thickness  also  increases  the  diffusion  resistance  in 
the  region  of  GDL  directly  in  contact  with  the  current  collector 
ribs.  Hence  there  is  an  optimum  thickness  below  which  the  cur¬ 
rent  density  reduces  for  smaller  GDL  thicknesses  [58].  In  our 
previous  work  [36],  we  showed  that  most  of  the  heat  generated 
is  transported  to  the  coolant  channel  through  the  current  collec¬ 
tor  ribs  and  a  very  small  fraction  of  the  total  heat  generated  is 
transported  to  the  gas  channel.  Smaller  GDL  thicknesses  reduce 
the  area  of  heat  transfer  from  the  high-temperature  region  (the 
region  of  the  GDL  directly  exposed  to  the  gas  channel)  to  the 
current  collector  ribs.  This  reduction  in  heat  transfer  area  results 
in  higher  temperatures  in  the  region  directly  exposed  to  the  gas 
channel  [36] .  In  the  region  of  the  GDL  directly  in  contact  with  the 
current  collector,  thinner  GDL  also  reduces  the  average  temper¬ 
ature  and  promotes  liquid  accumulation  in  this  region.  We  also 
showed  that,  using  GDLs  with  highly  anisotropic  thermal  con¬ 
ductivity  leads  more  uniform  temperatures.  Thus,  in  case  of  thin 
GDLs,  anisotropic  thermal  conductivity  plays  a  more  important 
role  in  maintaining  uniform  temperatures  and  reducing  the  risk 
of  liquid  accumulation  than  in  case  of  thicker  GDLs  [36],  We 
confirmed  the  conclusions  of  our  previous  work,  obtained  using 
a  single-phase  model,  using  the  present  two-phase  model. 

7.  Conclusions 

A  two-dimensional  two-phase  model  is  presented  to  pre¬ 
dict  the  effect  of  GDL  thermal  conductivity  on  the  polarization 
behavior  under  different  operating  conditions.  The  important 
conclusions  from  the  study  are  summarized  by: 

1 .  The  fuel  cell  generates  high  current  densities  with  moderately 
high  temperatures  for  a  low  value  of  through-plane  thermal 
conductivity  and  a  high  in-plane  thermal  conductivity  of  the 
GDL. 

2.  By  tailoring  the  GDL  to  have  high  through-plane  thermal 
conductivity  near  the  inlet,  membrane  dehydration  may  be 
avoided  for  low  humidity  operating  conditions. 

3 .  Similarly,  the  current  density  may  be  maximized  at  low  humi¬ 
dity  operating  conditions  by  tailoring  the  GDL  to  have  high 
through-plane  thermal  conductivity  near  the  inlet  and  pro¬ 
gressively  decreasing  through-plane  thermal  conductivity  at 
distances  farther  away  from  the  inlet  along  the  flow  channel. 

4.  Narrower  current  collector  ribs  lead  to  higher  current  densi¬ 
ties  for  fuel  cells  using  fully  humidified  cathode  inlet  flow. 

5.  Near  the  inlet,  fuel  cells  with  wider  ribs  produce  higher  cur¬ 
rent  density  than  those  with  narrower  ribs  if  the  cathode  inlet 
humidity  is  below  a  certain  value. 
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6.  The  current  density  can  be  maximized  at  lower  cathode  inlet 
humidities  by  using  wider  ribs  in  the  region  near  the  inlet  and 
narrower  ribs  in  the  region  away  form  the  inlet. 
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